Compute interaction between Sleep deprivation effect and genotype effect at QTL peaks

Read data

File<-"../Data/Metabolite_BXD.txt"
x<-read.table(File,header=T,sep="\t")
n<-x[1,] # Contains the metabolites type
names(n)<-gsub(".","_",names(n),fixed = TRUE) # reformate a bit the data
x<-x[-1,] # We remove the first lines containing the metabolites type

Some of our results contains Null values (-1), we replace them by NA values

head(x[1:10,1:10])
x[x==-1]<-NA

# Take only BXD
x<-x[-which(x$Group %in% c("B6D2F1","C57BL6","D2B6F1","DBA")),]

Read FC-mQTL results

fcmQTL<-read.table("../../QTL/Metabolite/FC/QTLlist.txt",header=T)

Get marker

markers<-read.table("Genotype.FormatedName.geno",header=T)
n<-markers$Locus
colnames(markers)
##  [1] "Chr"    "Locus"  "cM"     "Mb"     "BXD005" "BXD029" "BXD032"
##  [8] "BXD043" "BXD044" "BXD045" "BXD048" "BXD098" "BXD049" "BXD050"
## [15] "BXD051" "BXD055" "BXD056" "BXD061" "BXD063" "BXD064" "BXD065"
## [22] "BXD097" "BXD066" "BXD067" "BXD070" "BXD071" "BXD073" "BXD103"
## [29] "BXD075" "BXD079" "BXD081" "BXD083" "BXD084" "BXD085" "BXD087"
## [36] "BXD089" "BXD090" "BXD095" "BXD096" "BXD100" "BXD101" "B61"   
## [43] "DB1"
colnames(markers)<-gsub("BXD0","",colnames(markers))
colnames(markers)<-gsub("BXD","",colnames(markers))
colnames(markers)<-paste("BXD",colnames(markers),sep='')

markers<-markers[,colnames(markers) %in% unique(x$Group)]

# order marker
markers<-markers[,order(match(colnames(markers),as.character(unique(x$Group))))]

rownames(markers)<-n

Recompute quickly to get the peak

and then compute the interaction between genotype and sleep deprivation at QTL peak

library(qtl)
QTL<-read.cross("csv",file="../../QTL/Metabolite/FC/MetabolitesFC.QTLR",genotypes=c("BB", "BD", "DD","EE","FF"),allele=c("B", "D"),estimate.map=FALSE,na.strings=c("NA","-"))
##  --Read the following data:
##   33  individuals
##   11186  markers
##   126  phenotypes
## Warning in fixXgeno.f2(cross, alleles): --Omitting 75 male heterozygote
## genotypes on the X chromosome.
## Warning in summary.cross(cross): Some markers at the same position on chr
## 1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,X; use jittermap().
##  --Cross type: f2
QTL<-convert2risib(QTL)
## Warning in convert2risib(QTL): Omitting 6411 genotypes with code==2.
QTL<-jittermap(QTL,amount=1e-6)
QTL <- calc.genoprob(QTL, step=1,error.prob=0.05)

intpvals<-c()
effects<-c()
topmarkers<-c()

QTLinfo<-matrix(nrow=0,ncol=10)
SavedMarker<-matrix(ncol=33,nrow=0)
colnames(SavedMarker)<-colnames(markers)

markers<-markers[,order(match(colnames(markers),as.character(unique(x$Group))))]

for (k in as.character(unique(fcmQTL$Phenotype))){
  
  i<-which(fcmQTL$Phenotype==k & fcmQTL$Lodscore == max(fcmQTL[fcmQTL$Phenotype==k,"Lodscore"]))
  
  id<-fcmQTL[i,"PhenotypeNumber"]
  QTL.scan<-scanone(QTL,method="em",pheno.col=id)
  st<-which(rownames(QTL.scan) == fcmQTL[i,"MarkerStart"])
  ed<-which(rownames(QTL.scan) == fcmQTL[i,"MarkerEnd"])
  
  QTLscan_small <- QTL.scan[st:ed,]
  
  # Markers at QTL peak
  markersQTL<-QTLscan_small[QTLscan_small$lod==max(QTLscan_small$lod),]
  
  #remove pseudomarker
  gr<-grep("c\\d+\\.loc\\d+",rownames(markersQTL),perl=T)
  if (length(gr)>0){
    markersQTL_Nopseudo<-rownames(markersQTL)[-gr]
  }else{
    markersQTL_Nopseudo<-rownames(markersQTL)
  }
  
  #QTLscan_small<-QTLscan_small[- grep("c\\d+\\.loc\\d+",markersQTL,perl=T),]
  
  # if only pseudo marker available:
  if (length(markersQTL_Nopseudo) == 0){
    #estimate most probable genotype
    marker<-rownames(markersQTL)[1]
    
    Chrom<-fcmQTL[i,"QTL_Chromosome"]
    genoproDD<-pull.genoprob(QTL, chr=Chrom)[,paste(marker,':DD',sep='')]
    
    id<-as.character(QTL$pheno$Mouse_ID)
    gm<-rep("B",length(id))
    names(gm)<-id
    gm[genoproDD>0.5]<-'D'
    
  # A normal marker is available
  }else {
    marker <- markersQTL_Nopseudo[1]
    gm<-as.character(as.matrix(markers[marker,]))
    names(gm)<-colnames(markers)
  }
  
  if (nrow(QTLscan_small) > 0){
      maxid<-which(QTLscan_small[,"lod"]==max(QTLscan_small[,"lod"]))[1]
      marker<-rownames(QTLscan_small)[maxid]
      
      gms<-gm[order(match(names(gm),colnames(SavedMarker)))]
      SavedMarker<-rbind(SavedMarker,gms)
      rownames(SavedMarker)[nrow(SavedMarker)]<-marker

      geno<-gm[as.character(x$Group)]
    
      geno<-as.factor(geno)
      geno<-relevel(geno, ref="B")
      
      Pheno<-fcmQTL[i,"Phenotype"]
      Pheno<-gsub("_",'.',Pheno)
      idxmeta<-which(colnames(x) == Pheno)
      
      fit <-lm(as.numeric(as.character(x[,idxmeta]))~0+geno/factor(x$Group)+geno*factor(x$Condition))
      Bestim<-summary(fit)[[4]]["genoB","Estimate"]
      Destim<-summary(fit)[[4]]["genoD","Estimate"]
      SDestim<-summary(fit)[[4]]["factor(x$Condition)SD","Estimate"]
      DSDestim<-summary(fit)[[4]]["genoD:factor(x$Condition)SD","Estimate"]
      
      intpval<-summary(fit)[[4]]["genoD:factor(x$Condition)SD","Pr(>|t|)"]
      
      SD_B<-log2(Bestim+SDestim)-log2(Bestim)
      SD_D<-log2(Destim+SDestim+DSDestim)-log2(Destim)
      FCeff<-SD_B-SD_D
      
      effects<-c(effects,FCeff)
      intpvals<-c(intpvals,intpval)
      topmarkers<-c(topmarkers,marker)
  } else {
    intpvals<-c(intpvals,1)
    effects<-c(effects,0)
    topmarkers<-c(topmarkers,NA)
  }
  QTLinfo<-rbind(QTLinfo,fcmQTL[i,])

}
## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 1 individuals with missing phenotypes.
## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 1 individuals with missing phenotypes.
## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 2 individuals with missing phenotypes.
nrow(QTLinfo)
## [1] 89
QTL.scan[marker,]
intpvals<-p.adjust(intpvals,method="fdr")

QTLinfo$pvalint<-intpvals
QTLinfo$intlogFC<-effects
QTLinfo$topmarker<-topmarkers

fctest<-QTLinfo[,c(1,5,7,11,12,13)]
fctest[order(abs(fctest$intlogFC),decreasing = T),]
write.table(QTLinfo,file="Genotype_SD_Interaction_Metabolites.txt",sep="\t",quote=F,row.names=F)
plotFCMetabo<-function(WithVariant,WithoutVariant,Metabo,main){
  File<-"../../Rdata/Metabolites_BXD_alldata.txt"
  x<-read.table(File,header=T,sep="\t")
  x[x==-1]<-NA
  x<-x[-1,]
  m<-x[,-c(1,2,5)]
  m$Group
  #WithVariant<-c("DBA",WithVariant)
  #WithoutVariant<-c(WithoutVariant,"C57BL6")
  MetaNSDMean<-matrix(ncol=1,nrow=0)
  for (i in unique(m$Group)){
    MetaNSDMean<-rbind(MetaNSDMean,mean(as.numeric(as.character(m[m$Group==i & m$Condition=="NSD",Metabo])),na.rm=T))
  }
  rownames(MetaNSDMean)<-unique(m$Group)
  MetaSDMean<-matrix(ncol=1,nrow=0)
  for (i in unique(m$Group)){
    MetaSDMean<-rbind(MetaSDMean,mean(as.numeric(as.character(m[m$Group==i & m$Condition=="SD",Metabo])),na.rm=T))
  }
  rownames(MetaSDMean)<-unique(m$Group)
  FC<-MetaSDMean/MetaNSDMean
  d<-c()
  col<-c()
  geno<-c()
  for (i in rownames(FC)){
    if (i %in% WithoutVariant){
      d<-c(d,FC[i,])
      col<-c(col,rgb(0.5,0.5,0.5,0.8))
      geno<-c(geno,"B")
    }
    if (i %in% WithVariant){
      d<-c(d,FC[i,])
      col<-c(col,rgb(0.85,0.67,0.44,0.8))
      geno<-c(geno,"D")
    }
    if (i =="C57BL6"){
      d<-c(d,FC[i,])
      col<-c(col,rgb(0,0,0,1))
      geno<-c(geno,"B")
    }
    if (i =="DBA"){
      d<-c(d,FC[i,])
      col<-c(col,rgb(0.7,0.5,0.3,1))
      geno<-c(geno,"D")
    }
    if (i == 'D2B6F1' | i == 'B6D2F1'){
      d<-c(d,FC[i,])
      col<-c(col,rgb(1,1,1,1))
      geno<-c(geno,"BD")
    }
  }
  names(geno)<-names(d)
  
  barplot(height=log2(d[order(d)]),las=1,horiz=TRUE,col=col[order(d)],main=main,xlab="log Fold-Change")
}
QTLinfo<-QTLinfo[order(abs(QTLinfo$intlogFC),decreasing = T),]
par(mfrow=c(2,2))
for (i in 1:nrow(QTLinfo)){
  if (QTLinfo[i,"pvalint"]<=0.05){
      
      marker <- as.character(QTLinfo[i,"topmarker"])
      mdata<-SavedMarker[marker,]
      
      WithVariant<- names(mdata)[mdata=="D"]
      WithoutVariant<- names(mdata)[mdata=="B"]
      Metabo<-gsub('_','.',QTLinfo[i,"Phenotype"])
      main<-Metabo
      plotFCMetabo(WithVariant,WithoutVariant,Metabo,main)
    
  }
}